import numpy as np
import matplotlib.pyplot as plt


acelera_pendulo
euler_cromer
maxminv
medir_periodo
acc_toque
acc_i

#P10 - 1 
#a)movimento pêndulo método de Euler-Cromer

##b)a equação diferencial pode ser escrito -g/L * theta,
#expressao analítica é theta(t) = A * cos(sqrt(g/L) * t + phi)
#compare esta solução com a solução numérica obtido em alinea a)


#P10 - 2
#interpolação com polinómios de Lagrange período das oscilações

#P10 - 3
#a) período das oscilações
#b) gráfico do log do per 𝑇 em funçaõ do log do comp L
#c) ajuste linear, declive e erro


def acelera_pendulo(theta, L, g):
    # Aceleração do pêndulo
    a = -g / L * np.sin(theta)
    return a


def euler_cromer(theta, v, acelera, L, g, dt):
    # Calcula a aceleração
    a = acelera(theta, L, g)
    
    # Atualiza a velocidade
    v += a * dt
    
    # Atualiza o ângulo
    theta += v * dt
    
    return theta, v


theta_values = np.zeros(n)
v_values = np.zeros(n)
theta_values[0] = theta0
v_values[0] = v0

for i in range(1, n):
    theta_values[i], v_values[i] = euler_cromer(theta_values[i-1], v_values[i-1], acelera_pendulo, L, g, dt)



def maxminv(xm1,xm2,xm3,ym1,ym2,ym3):  
    # Máximo ou mínimo usando o polinómio de Lagrange
    # Dados (input): (x0,y0), (x1,y1) e (x2,y2) 
    # Resultados (output): xm, ymax 
    xab=xm1-xm2
    xac=xm1-xm3
    xbc=xm2-xm3

    a=ym1/(xab*xac)
    b=-ym2/(xab*xbc)
    c=ym3/(xac*xbc)

    xmla=(b+c)*xm1+(a+c)*xm2+(a+b)*xm3
    xm=0.5*xmla/(a+b+c)

    xta=xm-xm1
    xtb=xm-xm2
    xtc=xm-xm3

    ymax=a*xtb*xtc+b*xta*xtc+c*xta*xtb
    return xm, ymax




# Use interpolação com polinómios de Lagrange para medir o período das oscilações

def medir_periodo(theta, t):
    maximos_tempos = []
    for i in range(1, len(theta)-1):
        if theta[i-1] < theta[i] and theta[i] > theta[i+1]:
            t_max, _ = maxminv(t[i-1], t[i], t[i+1], theta[i-1], theta[i], theta[i+1])
            maximos_tempos.append(t_max)
    if len(maximos_tempos) >= 2:
        return maximos_tempos[1] - maximos_tempos[0]
    else:
        return None



# Verifica se encontrou pelo menos dois máximos
if len(maximos_tempos) >= 2:
    T_medido = maximos_tempos[1] - maximos_tempos[0]
    T_teorico = 2 * np.pi * np.sqrt(L / g)

    print(f"Período medido: {T_medido:.4f} s")
    print(f"Período teórico: {T_teorico:.4f} s")
    print(f"Erro relativo: {abs(T_medido - T_teorico)/T_teorico * 100:.2f}%")
else:
    print("Não foram encontrados dois máximos para calcular o período.")



def acc_toque(dx,d):
    # calcular a aceleração de uma esfera devido ao contacto com a esfera à sua direita
    k = 1e7
    q = 2.0
    if dx<d:
        a = (-k*abs(dx-d)**q)/m
    else:
        a = 0.0
    return a


def acc_i(i,x):
    #calcular a aceleração de esfera i, cuja posicao de equilibrio é d*i
    a = 0

    if i>0: # a primeira esfera não tem vizinho à sua esquerda
        a -= acc_toque(x[i] - x[i-1],d)
    if i < (N-1): # a última esfera não tem vizinho à sua direita
        a += acc_toque(x[i+1] - x[i], d)

    # aceleração de gravidade, afeta todas as esferas
    a -= g*(x[i]-d*i)/l
    return a